Collapse instability of solitons in the nonpolynomial 
Schrodinger equation with dipole-dipole interactions 

G. Gligoric^ A. Maluckov^, Lj. Hadzievski^ , and B. A. 
Malomed^ 

Abstract. A model of the Bose-Einstein condensate (BEC) of dipolar atoms, 
confined in a combination of a cigar-shaped trap and optical lattice acting in the axial 
direction, is studied in the framework of the one-dimensional (ID) nonpolynomial 
Schrodinger equation (NPSE) with additional terms describing long-range dipole- 
dipole (DD) interactions. The NPSE makes it possible to describe the collapse of 
localized modes, which was experimentally observed in the self-attractive BEC confined 
in tight traps, in the framework of the ID description. We study the influence of the DD 
interactions on the dynamics of bright solitons, especially as concerns their collapse- 
induced instability. Both attractive and repulsive contact and DD interactions are 
considered. The results are summarized in the form of stability/collapse diagrams 
in a respective parametric space. In particular, it is shown that the attractive DD 
interactions may prevent the collapse instability in the condensate with attractive 
contact interactions. 
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1. Introduction 

In the mean-field approximation, the dynamics of Bose-Einstein condensates (BECs) 
obeys the 3D Gross-Pitaevskii equation (GPE) [1], from which an effective ID equation 
can be derived, in various settings, for condensates trapped in prolate traps [2]-[7j. In the 
simplest case, which corresponds to a sufficiently low BEG density, the reduction of the 
3D equation for the BEG trapped in the "cigar-shaped" configuration leads to the one- 
dimensional (ID) cubic nonlinear Schrodinger equation (NLSE) [8]. The main restriction 
on the use of this equation is its inability to describe the onset of collapse of localized 
states, which was theoretically predicted in the 3D setting and experimentally observed 
in the self-attractive BEG [9], [10]. However, without imposing the constraint of a very 
low density, the reduction of the 3D GPE leads to the ID equation with a nonpolynomial 
nonlinearity, alias the nonpolynomial Schrodinger equation (NPSE) [U [TT]. The ID 
NPSE with the attractive sign of the nonlinearity enables the description of the collapse 
dynamics and produces results which are corroborated by direct simulations of the 
underlying 3D GPE [12]. 

On the other hand, it is known that the mean-field description of BEGs trapped in 
a very deep optical-lattice (OL) potential can be well described by the corresponding 
discrete equations. In particular, discrete forms of the ID GPE with the cubic 
nonlinearity [131 El US], and of the ID NPSE [16] have been studied in detail. Basic 
features of the ID continual equations describing BEGs trapped in deep OLs find their 
counterparts in the discrete models. In particular, the ability of the continual ID 
NPSE to capture the onset of the collapse is also shared by the corresponding discrete 
equations. 

A new variety of the BEG dynamics, which is dominated by long-range (nonlocal) 
interactions, occurs in dipolar condensates, which may be formed by magnetically 
polarized ^^Gr atoms [17], dipolar molecules p^, or atoms with electric dipole moments 
induced by an external field [19]. A review of dynamical effects produced by the dipole- 
dipole (DD) interactions in condensates can be found in Ref. [20]. In particular, 
conditions for the stability of the trapped dipolar BEG against collapse were studied 
in detail [21]. A possibility of the creation of 2D solitons in dipolar condensates was 
predicted too. Namely, isotropic solitons [22] and solitary vortices [23] may exist if the 
sign of the dipole-dipole (DD) interaction is inverted by means of rapid rotation of the 
dipoles [21]. On the other hand, stable anisotropic solitons can be supported by the 
ordinary DD interaction, if the dipoles are polarized in the 2D plane [25]. Solitons 
supported by nonlocal interactions were also predicted and realized in optics, making 
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use of the thermal nonhnearity [20] • 

A natural extension of the consideration of the dipolar BEC includes the OL 
potential, which, in the discrete limit, leads to the model with the long-range DD 
interactions between lattice sites [271 [28| [29] . Very recently, ID solitons have been 
studied in the framework of the continual cubic GPE, assuming the competition of local 
and nonlocal DD interactions, with or without the OL potential [30]. In the discrete 
limit, ID solitons supported by the DD interactions were recently studied too, in models 
with both the cubic [27] and nonpolynomial [28] onsite nonhnearity. In particular, the 
latter work predicts a possibility of suppressing the collapse by means of the long-range 
DD forces. The main purpose of the present work is to investigate the influence of the 
DD interactions on the collapse dynamics in the ID continual NPSE, which is the most 
adequate setting for the study of the onset and suppression of the collapse in the dipolar 
BEC loaded into a cigar-shaped trap. 

The paper is structured as follows. The model equation which includes the 
nonpolynomial local nonhnearity and nonlocal DD interactions, is formulated in Section 
II, where we also outline numerical techniques that we use in this work. Focusing on 
the study of fundamental bright solitons, we report basic results for their existence and 
stability in Section III. The core part of paper is Section IV, where we study the influence 
of the DD interaction on the solitons' collapse. The paper is concluded by Section V. 



2. The model 



The dynamics of a BEC at zero temperature is accurately described by the 3D Gross- 
Pitaevskii equation (3D GPE) [Ij. When the condensate is confined by a harmonic 

1 /2 

potential with frequency u± and respective length a± = {h/mu±) in the transverse 
plane, and by generic potential V (z) in the axial direction, it was shown in Ref. [1] 
that the corresponding 3D GPE can be reduced to the ID NPSE for wave function 
ip{z,t), which is subject to normalization \ip {z)f dz = 1. The equation includes 
the OL potential, with depth Vq and wavenumber K, and the long-range DD (cf. the 
ID equations introduced in Refs. [271 130] ): 



Here, 7 = —INa^ muj^j h is the effective strength of the local interaction, with N 
the total number of atoms in the condensate, and the scattering length of atomic 
collisions (a^ < corresponds to attraction) [4j. Further, G = g {1 — Scos^O) is the 
coefficient which defines the DD interaction, where g is a. positive coefficient, and 6 the 
angle between the z axis and the orientation of the dipoles. 

Replacing wave function ip hy f = ^/ItIV') transform Eq.([T]) into a normalized 
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form, 

at zaz ^l_^|j|2 

where K = sgn(7) is the sign of the local interaction (K = +1 corresponds to the 
attraction), and V = G/ I7I measures the relative strength of the DD and contact 
interactions. 

One case of obvious interest for the ID setting is that when the dipoles are aligned 
with the z axis, i.e., ^ = and F = —2g/ \'~f\ < 0, which means that the DD interaction 
is attractive. Another relevant case corresponds to the dipoles oriented perpendicular to 
the z axis, i.e. 6 = tt/2 and T = g/ > 0, which implies the repulsive DD interaction. 
Thus, the sign of F in the present model defines the character of the DD interaction. 

Stationary solutions to Eq. ([2]), with chemical potential /i, are sought for as 
/ = uexp{—ifit), with function u satisfying the stationary equation, 

l^u= - 9 ^ + n) sm [Kz) u + --U 

\j\-^\uf 

+ Fn(z) M^dz'. (3) 
J —00 \z I 

To present the results for soliton families, we will use the norm defined by 

/+00 
\u{z)fdz (4) 

(according to the above definitions, P is identical to I7I , but parameter 7 is not used 
below) . 

Experimentally adjustable coefficients in this model are the relative strength 
of the DD/contact interactions, F, and the norm of the stationary wave function, 
P ^ N \ as \ ^/muj^jli [m [ini [27]. In particular, P ~ 1 corresponds to ~ 1000 atoms in 
the ^^Cr condensate [HI [16]. Actually, F can be made both positive and negative, and its 
absolute value may be altered within broad limits by means of the Feshbach resonance 
(without the application of the Feshbach resonance, |F| ~ 0.15 in the condensate of 
chromium atoms) [171 [31] . 

Stationary equation (j3]) was solved numerically by means of an algorithm based on 
the shooting method. We restrict the analysis to fundamental solitons with a single 
maximum. Therefore, we chose the following boundary conditions in the shooting 
procedure: the first derivative of the wave function must be zero at the point of the 
maximum amplitude, and the wave function must decay exponentially at infinity. Note 
that the integral term in Eq. ([3]) has the form of a convolution [32] , 
30 

V{z - z')\u{z')f dz! = V{z) X \u{z)f , (5) 
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hence its Fourier transform can be calculated as a product of the Fourier images of the 
interaction potential and local density, To calculate the integral term more 

accurately, we adopted an iterative shooting procedure. For the first iteration, we 
assumed that the potential of the interaction in the Fourier space is approximated by 
a constant [32] ■ In each next iteration, we used the previous solution to calculate the 
convolution integral numerically ([5]). 

Time- dependent equation ([2]) was solved by means of the split-step Fourier method 
[33]. To that end, Eq. ([2]) was cast into an operator form, 
df .... 

-^ = i{V + D + L + N)f, (6) 
where V = -Vq sin^ (Kz), D = -T \f{z')\'^/\z - z'\^dz\ L = dy{2dz'^), and N = 

— 1 /2 

- (1 - K Iff) ' [1 - (3/2)N Iff] . The split-step Fourier method was implemented 
through independent actions of all of the operators, V, D, L and N, within time step r 
of the integration scheme, f(t + T, z) = e~*'^'^e~*'^^e~*'^^e~*'^^/(t, z) (we fixed r = 10""^). 

3. Fundamental solitons 

Stationary wave function u{z) was found from the solution of stationary equation ([3]). 
In the case of the noninteracting condensate (7 = 0), and without the DD interactions 
(r = 0), equation ([3]) is tantamount to the Mathieu equation, 

d'^u 

— -[2qcos{2y)-p]u = 0, (7) 

where y = Kz, q = —Vq/2K^, and p = 2 {q — 1/ K"^ + fi/K"^), hence the well-known 
bandgap diagram for the Mathieu equation (see, e.g., Ref. |33J) can be used. That 
diagram, translated into the notation adopted above, is displayed in Fig. [TJ The 
nonlinear localization of matter waves in the form of gap solitons may occur in gaps 
of the linear spectrum. In the case of the local attractive nonlinearity, fundamental 
solitons populate the semi-infinite gap. 

Figure El displays typical examples of fundamental solitons, with equal values of the 
norm, which were obtained as numerical solutions to Eq. ([3]) in cases of zero, attractive, 
and repulsive DD interactions. In the presence of the attractive DD interaction, the 
solitons are narrower and feature a higher amplitude, while in the case of the DD 
repulsion, they are wider, and have a lower amplitude than the respective soliton in the 
absence of the DD interaction. The peculiarity of the NPSE in the case of the attractive 
contact interaction is that the amplitude of the soliton is limited by a critical value. 
The collapse setting in when the amplitude attains that value [H [12]. In contrast to 
that, in the ordinary ID GPE the amplitude may grow indefinitely without causing the 
collapse. 

If the local interaction is repulsive (N = —1), fundamental solitons in the NPSE 
may only be supported by the attractive DD interaction, which must dominate over 
the contact repulsion. In that case, the difference between the NPSE and ordinary 
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ID GPE is not significant. In particular, as the collapse instability is absent, unstable 
fundamental solitons can only evolve into breathers. The situation in the case of the 
repulsion is similar to that in the case of the discrete NPSE |27j . 

4. The influence of the dipole-dipole interactions on the collapse 

In this section we focus on the dynamics of fundamental solitons in the presence of the 
attractive contact interaction (K = +1) and nonlocal DD of either sign. The collapse 
occurs in this case, as well as in the respective discrete setting [27] . 

A global characteristic of soliton families is the -P(/i) dependence, the norm versus 
the chemical potential. In Fig. [3] we display the -P(/i) curves obtained for different values 
of r. The stability of the soliton may be assessed according to two different conditions: 
dP/d^ < [the Vakhitov-Kolokolov (VK) criterion], and the absence of eigenvalues for 
small perturbations with positive real parts (the spectral condition), see Ref. [35] and 
references therein. 

In the presence of the attractive contact interaction, the P{fi) dependences feature 
two different regions, one where the VK criterion is fulfilled, and another one, where it is 
violated. The region where the VK criterion holds expands/shrinks with the increase of 
the strength of the attractive/repulsive DD interactions. These results are in qualitative 
agreement with findings for on-site solitons in the discrete version of ID NPSE [27]. 

The spectral stability condition was examined by numerically computing the 
corresponding eigenvalues, using linearized equations for small perturbations. It has 
been found that the spectral condition is violated in the entire parameter space which 
was explored, see Fig. HI The eigenvalue spectra are quite similar to those found for 
inter-site solitons in the discrete NPSE with the DD interaction, at large values of the 
inter-site coupling constant C, which corresponds to the quasi-continuum limit [27j . 
Note that the spectra for solitons in the NPSE exhibit an abrupt (quasi-exponential) 
growth of real parts of the eigenvalues [see Fig. |4](a)], which does not happen in the 
respective GPE with the cubic nonlinearity, cf. Fig. Hj^b). An interesting fact revealed 
by the numerical analysis is that the threshold of the abrupt growth of the real part 
of the complex eigenvalues corresponds to the onset of the collapse instability of the 
solitons. Indeed, direct simulations of Eq. ([2]) show that unstable solitons with /i 
higher than the threshold value evolve into robust localized breathers [Fig. [5]^a)], while 
the solitons with /i taken below the threshold exhibit the collapse instability, which 
manifests itself through simultaneous decrease of the soliton's width and growth of the 
amplitude towards the limit value, as seen in Fig. [5t^b). Points of the onset of the 
collapse instability are also marked in Fig. [3l It is worthy to note that the instability 
sets in earlier than the norm attains the maximum values (full curves P{n) could be 
drawn in spite of the instability, as they were obtained from the numerical solution of 
stationary equation ([3])). 

With the increase of the strength of the attractive DD interaction, the threshold 
value of n becomes lower, hence the region in the parameter space where unstable 
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solitons do not collapse but rather evolve into breathers expands. These results are 
presented as the collapse diagram in the (r,/i) plane (see Fig. [6]), where the border line 
separates the collapse region from that where the formation of breathers takes place. 
Qualitatively, the nearly linear dependence of the threshold on T{n) can be understood, 
assuming that the nonlocality range is very large, covering the entire soliton. Indeed, 
in that case the last term in Eq. ([3]) is proportional to TPu{z), which implies a linear 
shift of /i. 

Thus, the collapse instability of the fundamental solitons may be suppressed by 
using sufficiently strong attractive DD interactions. This prediction of the stability 
analysis based on the computation of the eigenvalues was fully corroborated by direct 
simulations of Eq. ([2]) . 

5. Conclusion 

The purpose of this work was to achieve a better understanding on the influence of 
the nonlocal DD (dipole-dipole) interactions on the stability and collapse of localized 
nonlinear modes in the BEC trapped in a combination of a tight transverse parabolic 
potential and a relatively loose periodic OL potential acting in the axial direction. 
To this end, we have introduced the model based on the one-dimensional NPSE 
(nonpolynomial Schrodinger equation), which includes the contact and DD interaction 
terms, as well as the OL potential. Both attractive and repulsive signs of the contact 
and DD interactions were considered. The analysis was focused on fundamental solitons 
in the semi-infinite gap. While all the stationary solitons are unstable, an essential 
conclusion is that the attractive DD interactions may prevent the collapse instability of 
the fundamental solitons, replacing the onset of the collapse by the transformation of 
the solitons into robust breathers. This general result is consistent with findings recently 
reported for the ID discrete version of the NPSE. 
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